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ABSTRACT 

We investigate the molecular evolution and D/H abundance ratios that de- 
velop as star formation proceeds from a dense-cloud core to a protostellar core, 
by solving a gas-grain reaction network applied to a 1-D radiative hydrodynamic 
model with infalling fluid parcels. Spatial distributions of gas and ice-mantle 
species are calculated at the first-core stage, and at times after the birth of 
a protostar. Gas-phase methanol and methane are more abundant than CO 
at radii r < 100 AU in the first-core stage, but gradually decrease with time, 
while abundances of larger organic species increase. The warm-up phase, when 
complex organic molecules are efficiently formed, is longer-lived for those fluid 
parcels in-falling at later stages. The formation of unsaturated carbon chains 
(warm carbon-chain chemistry) is also more effective in later stages; C + , which 
reacts with CH4 to form carbon chains, increases in abundance as the envelope 
density decreases. The large organic molecules and carbon chains are strongly 
deuterated, mainly due to high D/H ratios in the parent molecules, determined 
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in the cold phase. We also extend our model to simulate simply the chemistry in 
circumstellar disks, by suspending the 1-D infall of a fluid parcel at constant disk 
radii. The species CH3OCH3 and HCOOCH3 increase in abundance in 10 4 — 10 5 
yr at the fixed warm temperature; both also have high D/H ratios. 

Subject headings: stars: formation, ISM: clouds, ISM: abundances 



Introduction 



In star-forming cores, temperature and density vary both temporally and spatially. 
The rates of chemical processes, which determine the molecular abundances in these cores, 
also change according to the physical conditions. In prestellar cores with high densities 
and low temperatures, various atoms and molecules are depleted onto grains, where they 
are hydrogenated by grain-surface reactions. In the later, protostellar phase, on the other 
hand, molecules are desorbed from the grains back into the gas phase in the central high- 
temperature regions. The molecular abundances in the cores thus change dramatically in 
the course of star formation. 

It is well established that molecular clouds are generally not in chemical equilibrium; the 
chemical timescale is comparable or longer than the dynamical timescale of the clouds. For 
example, the adsorption timescale of gaseous species onto grain surfaces in cold prestellar 



cores is 
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where a and n grain are the radius and number d ensity of the grain particles, respectively, and 



S is the sticking probability (e.g. lHerbstlll993l ). At the typical density of molecular clouds, 
% ~ 10 4 cm " 3 , this timescale is comparable to the free-fall timescale of a cold core (e.g. 



Spitzerill978f ) 
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Thus, hydrodynamic-chemical models are required to investigate theoretically the abun- 
dances and spatial distributions of molecules in star-forming cores. 

Since the hydro-chemical models combine hydrodynamics and chemical reaction-network 
models, there is a variety of previous work related to this topic. The chemistry in prestel- 
lar cores ha s been investigated using Bonner-Ebert s pheres and/or the isothermal collapse 



model (e.g. lAikawa et al.l l200ll ; iKeto fc CaseUil 120101 ). The chemistry in protostellar cores 
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is often investigated using a pseudo-time-dependent gas-phase chemistry model with a con- 
stant temperature (~ 200 K) and density (rin ~ 10 6 cm" 3 ), in which a high abundance of 



ice sublimates is given as an initial condition (e.g. iCharnley et al.lll992f ). While such simple 
models give chemical timescales of various species, which is useful in the interpretation of 
protostellar cores, ice abundances should ideally be determined by the gas-grain chemical 
model in the prestellar phase. 



Rodgers fc Charnlev (|2003l) inv estigated the chemistry in a protostellar core using the 



inside-out collapse model ( jShu Ill977[ ); they calculated radial distributions of molecular abun- 
dances at 10 2 , 10 3 , 10 4 and 10 5 yr after the collapse of the singular isothermal sphere. Tempo- 
ral variation of the temperature di s tribut ion, which is a critical parameter for the chemistry, 
was adopted from lAdams & Shu I (119851 ). The model shows the onion-skin type structure, 
in which more volatile species are sublimated to the gas phase at outer radii, with the sub- 
limation radius moving outwards as the core becomes hotter. The model did not include 
a prestellar phase, since the inside-out collapse model corresponds to the main accretion 
phase, and the ice abundance was set as an initial condition. 



Lee et al.l ( 120041 ) went one step farther; they solved the chemistry in a core starting from 



a prestellar core and evolving to a protostellar core self-consistently. Their physical model 
of the core is a combination of a series of Bonner-Ebert spheres and inside-out collapse. 
The spatial and temporal variation of the temperature in the protostar is calculated by 
radiation transfer, assuming the luminosity of the central star as a function of time. The 
model clearly shows that the sublimation of ice (e.g. CO) affects the gas-phase chemistry 
significantly. Their chemical model, however, does not include grain-surface reactions except 
for H2 formation; dust-grain ice mantles are composed of molecules accreted directly from 
the gas phase. 



Garrod fc Herbstl (120061 ). on the other hand, constructed a comprehensive gas-grain 
chemical network with various grain-surface reactions. Their chemical network is solved 
with simple physical models in which the density is constant and the temperature increases 
with time. They showed that some large organic species are efficiently formed on grain 
surfaces starting at temperatures of 30-40 K and sublimated at higher temperatures. Their 
model accounts for the high abundance of large organic molecules such as methyl formate 
(HCOO CH3) and dimethyl et he r (CH3OCH3) in low -mass protostellar cores, known as hot 



corinos (Ceccarelli et al. 2007). Aikawa et a 



chemical network of Garrod Sz Herbst 



(120081 ). hereinafter Paper I, then applied the 
(2006) to a 1-D radiation- hydrodynamic model of a 



star-forming core by iMasunaga fc Inutsukal (120001 ) . As expected, organic species such as 
CH3CN and HCOOH are formed in the warm-up phase and become abundant in the central 
region (T > 100 K) of the protostellar core. In addition, the model shows that (unsaturated) 
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carbon-chain species increase in abundance inwards at T ~ 25 K, which is observed as the 
so-called "Warm Carbon Chain Chemistry (WCCC)" JSakai et al.llioosh . 

In this work, we update and improve the model of Paper I. The aim of the work is to 
investigate 



1. the spatial distribution of molecules at several evolutionary stages from a prestellar 
core to a protostellar core, 

2. the evolution of molecular D/H ratios, 

3. how the chemical composition in a disk differs from that in cores. 

Paper I showed spatial distributions of molecules only at 9.3 x 10 4 yr after the protostellar 
birth, although the model starts from a prestellar core. Here we show molecular distributions 
at assorted evolutionary stages to see how the distributions and peak abundances of hot 
corino species and carbon chains vary with time. We also extend the chemical reaction 
network to include both singly- and multiply-deuterated species. It is well established that 
molecular D/H ratios are enhanced by exothermic exchange reactions at low temperature. 
Once the star is formed, the core becomes warmer and the backward reactions (endothermic- 
exchange reactions) become efficient. How the D/H ratios vary after protostellar birth is of 
importance to link the observed high D/H ratios in the ISM to the isotopic ratios observed 
in planetary matter such as comets and meteorites. In order to investigate such a link, it 
is also of great interest to determine how the chemical composition in a disk differs from 
that of cores. For the material in an infalling envelope, the duration of the warm-up phase 
and warm chemistry is determined by the size of the warm region divided by the free-fall 
velocity, which is rather short, as shown in Paper I. In circumstellar disks, on the other 
hand, material can stay in warm regions for longer timescales, since the disk is supported 
by rotation, at least partially. Because our model is spherical and thus does not include 
the structure and dynamics of disk formation, we briefly investigate the effect of prolonged 
warm temperature chemistry in the disk by calculating pseudo-time-dependent models with 
the initial abundances set by the infalling core model. 

In the following section, we briefly describe our physical model of the star-forming core 
and the chemical reaction network (§2). In §3, we report the spatial distribution of molecular 
abundances and D/H ratios at assorted evolutionary stages. In §4, we discuss the effect on 
the chemistry of 2D structure in the core; we calculate a pseudo-time-dependent model of 
a warm disk-like phase. We also compare our model results with observational studies, and 
discuss the effect of grain-surface abstraction and substitution reactions on D/H ratio of 
methanol. We summarize our conclusions in 85. 
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Model 



2.1. Physical model of a star-forming core 



As in Paper I, we adopt the ID (spherical) radiatio n hy drodynamic model of low-mas s 
star formation by lMasunaga. Miyama fc Inutsukal (119981 ) and iMasunaga &: Inutsuka] ( 120001 ) . 
Initially, the central density of the molecular cloud core is p = 1.4 x 10 -19 g cm -3 , which 
corresponds to a number density of hydrogen nuclei % ~ 6x 10 cm -3 . The outer boundary 
is fixed at r = 4 x 10 4 AU, so that the total mass is 3.852 M , which exceeds the criti- 
cal mass for gravitational instability. The contraction is almost isothermal as long as the 
cooling rate overwhelms the compressional heating, but eventually the latter dominates and 
the temperature rises in the central region. Afterwards, the first core, which is a hydro- 
static core of H 2 gas, is formed. When the central temperature reaches ~ 2000 K (several 
hundred years after the first core formation), H 2 starts to dissociate and the first core col- 
lapses. The central region becomes hydrostatic again when the dissociation and ionization 
of hydrogen are completed; this signifies the birth of a protostar. In our model, the prestel- 
lar core evolves to the protostellar core in 2.5 x 10 5 yr. After the birth of the protostar, 
the model further follows the evolution for 9.3 x 10 4 yr, during which the protostar grows 
by mass accretion from the envelope. At each evolutionary stage, the model gives the to- 
tal luminosity of the core and the radial distribution of density, t emperature, and infall 
veloci ty at r > 10 -4 AU self-consistently. The original model by IMasunaga fc Inutsuka 



(200 0j) included wavelength-dependent radiation trans fer involvin g the opacities of the gas 



Iglesias fc Rogers Ill996t Alexander fc Ferguson 1 1 19941 ) and dust (IPreibisch et al. Ill993[ ) to 



obtain both the temperatures of the gas and the various components of dust (silicate, carbon 
and ice) separately. For simplicity, we adopt their gas temperature and assume that it is 
equal to the dust temperature, since the differences between the gas and dust temperatures 
are not siginificant in the region we are interested in. Figure Q] (a — c) shows the distribu- 
tion of density, temperature, and infall velocity at assorted evolutionary stages, in which we 
present the molecular distribution in §3: t = —5.6 x 10 2 yr, 4.3 x 10 2 yr and 9.3 x 10 4 yr. 
Here we define t = as the moment of the birth of a protostar. The total luminosity of the 
core is about 5.9 x 10~ 3 L Q at t = —5.6 x 10 2 yr, 0.88L© at t = 4.3 x 10 2 yr, and 24L Q at 
9.3 x 10 4 yr. The first core is formed at around t = —5.6 x 10 2 yr; Figure [T] (a — c) shows 
that at this time a hydrostatic core of a few AU size is formed. It should be noted that the 
density in the envelope decreases with time after the protostar is formed, while the density 
increases w ith time in the prestellar stage. More d e tailed explanations of core ev olution can 
be found in IMasunaga. Miyama &: Inutsukal (119981 ). IMasunaga &: Inutsukal (120001 ) and Paper 
I. 



Figured] (d — e) shows the temporal variation of density and temperature in fluid parcels 
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that fall tor = 2.5 AU at t = -5.6 x 10 2 yr, 4.3 x 10 2 yr and 9.3 x 10 4 yr. Naturally, the 
density and temperature increase with time. As can be seen in Eq. (2), the free-fall timescale 
is proportional to p~ x l 2 . While the core temperature rises in the central region, the infall 
velocity increases toward the center until the fluid parcel hits the surface of the central 
hydrostatic core (the first core or the second core). Hence the temporal variation of density 
and temperature accelerates in the infalling fluid parcels. In order to highlight the rapid rise 
near and at the final stage, the horizontal axis in Figure Q] (d — e) is set to be the logarithm 
of U 5Au — t, where t2.5Au corresponds to t = —5.6 x 10 2 yr for the dotted line, 4.3 x 10 2 yr 
for the dashed line, and 9.3 x 10 4 yr for the solid line. 

The temporal variation of temperature is critically important for chemistry. Garrod 
& Herbst (2006) show that the grain-surface reactions of heavy-element radical species are 
efficient in a temperature range of ~ 30 — 100 K. Once the temperature gets higher than 
~ 100 K, most volatile species are sublimated to the gas phase, where they begin to undergo 
gas-phase reactions. In Figure [1] (d — e), we can easily read the timescale of these warm 
phases. For example, the fluid parcel that reaches r = 2.5 AU at t = 9.3 x 10 4 yr stays at a 
temperature of 30 — 100 K for ~ 10 4 yr. After the fluid parcel gets warmer than 100 K, it 
falls to the central star in a few 100 yr. The fluid parcels that reach r = 2.5 AU at earlier 
time spend even less time in the warm regions. 



In the original model by iMasunaga fc Inutsukal (120001 ). the core starts contraction im- 



mediately, since the core mass exceeds the critical mass for gravitational instability. For 
our chemistry model, however, we need to set up the initial molecular abundances of this 
dense molecular core by assuming a simple history of fluid parcels before collapse. Since the 
formation of molecular cloud cores is out of the scope of this work, we simply assume that 
the core keeps its hydrostatic structure for 1 x 10 6 yr, implicitly assuming that turbulence 
supports it. In addition, we assume that the temperature in the pre-collapse phase is 10 
K, as typically observed in molecular clou ds. After 1 x 10 6 yr, we low er the temperature 



smoothly but rapidly to the initial value in IMasunaga fc Inutsukal (120001 ) , ~ 7 K 



2.2. Chemical reaction network 



We solve the rate equations of the gas-grain reaction network of iGarrod fc Herbstl (120061 ) 
under the time-dependent physical conditions of each fluid parcel, to obtain the radial distri- 
butions of molecular abundances both in the gas phase and ice mantle s at assorted evolution- 
ary stages. We use the gas-grain code Nautilus ( IHersant et al. 1 12009c ISemenov et al. 1120101 ) 
to solve the rate equations. A layered structure of the ice mantle is not considered (see §3.2). 
The chemical reaction network model and parameters are essentially the same as in paper 
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I, except for the following three updates: (i) Some gas-phase reaction rates are updated 



based on the KIDA database (http://kida.obs.u-bordeauxl.fr/). (ii) We included reaction s 
between carbon-ch ain species with H atoms or H 2 molecules, based on lHarada et al.l ( j2010l ). 
Hassel et al.l (j201l[ ) showed that these reactions significantly change the evolution of radical 
carbon-chain species at temperatures > 100 K. We do not, howeve r , disc riminate between 
isomers such as C-C3H2 and H 2 CCC, which was done by lHassel et al.l ( 1201 ll ). (iii) We include 
multi-deuterated species according to the following procedures. 

First, the species list of the network model is extended to include mono-, doubly-, and 
triply-deuterated species. Then the reaction network, both in the gas phase and on grain 
surfaces, is extended to include these deuterated species. For example, if we have the reac- 
tion AH + + BH — 7- A + BH^ in the original reaction list, we add 
AD+ + BH -> A + BHD+ 
AH+ + BD ->• A + BHD+ 
AD+ + BD ->• A + BD+ 

assuming that the rate coefficients are the same as for the original reaction. If there is more 
than one set of possible products, 
AHD + BD ABH + D 2 
AHD + BD ABD + HD, 

for example, we set the branching ratio statistically (e.g. 



Furuya et al.ll201ll ). In the above 



example, the branching ratio for the former and latter product channels is 1:2. We also 
assumed statistical branching ratios in dissociative recombination reactions, except for the 
recombination of deuterated H^; this assumption can result in underestimates of the molec- 
ular D/H ratio by a factor of two or so, since channels with a hydrogen atom as a product 
are more likely than ana logous ones with a D atom due to the greater speed of the H atom 
( iLe Petit k, Rouefil 120031) . For the recombinatio n of deuterated H^, we adopted branching 
ratios listed in iRoberts. Herbst fc Millar! (120041 ) . Indeed, these branching ratios are higher 
for reactions producing an H atom than the statistical value. For simplicity, we do not follow 
the position of the deuterium atoms in the reactions which include deuterated isomers. For 
example, in the reaction of H2CO + OD, the branching ratio of the three product channel 
H 2 CO + OD HCOOH + D 
H 2 CO + OD DCOOH + H 
H 2 CO + OD HCOOD + H 

is 1:1:1 in our current model, although the last branch would be the most probable in reality. 
So we do not aim to predict the abundance ratios of deuterated isomers. 



Finally we included reactions listed in lMillar. Bennet fc Herbstl (119891 ) and lRoberts. Herbst fc Millar 



( 120041 ): these are mostly exchange reactions that cause the initial isotopic fractionation. At 
low temperatures, of several tens of Kelvin or less, deuterium fractionation is mainly trig- 
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gered by Hjj~ + HD — > H2D + + H2 (and its multi-deuterated analogues), which is exothermic 
by 230 K. The effective rate of the backward reaction depends on the ortho/para ratio of 
H 2 , since the internal energy of ortho H 2 helps o yercome the endothe rmicity. But the or- 
tho/para ratio of H 2 is a complicated issue (e.g. iFlower et al. II2006I ). and there is so far 



no observational e viden ce for a high abundance of ortho H 2 in molecular clouds. Recently, 
Watanabe et al. I ( 120101 ) measured the spin temperature of H 2 formed on ice in the labora- 
tory, obtaining a value greater than approximately 200 K. But they also found that ortho H 2 
that is trapped on the ice surface is converted to para H 2 after formation. Such conversion 
of ortho to para has never been considered in the prediction of the o/p ratio of H 2 . In this 
work, we do not consider the o/p ratio; i.e. we assume the backward reaction of + HD 
is endothermic by 230 K. 

Inclusion of multi-deuterated species naturally increases the number of species and re- 
actions. In total, our model includes 1564 species and 45371 reactions, while the original 
model includes 655 species and 6309 reactions. 



As elemental abundances, we adopt the so-called "low-metal" values (Table 1 of lAikawa et al. 



( 120011 )). The species are assumed to be initially in the form of atoms or atomic ions except 
for hydroge n, which is entirely i n its molecular form. The elemental D/H ratio is set to 
1.5 x 10~ 5 (IPiskunov et al.lll997l ). All deuterium resides initially in the HD molecule. The 



cosmic-ray ionization rate for H 2 is set to 1.3 x 10 s . 

The sticking probability for gaseous neutral species onto grains is set to 0.5, indepen- 
dent of temperature. Although the probability may be lower at high temperatures (T > 
100 K), the value of the sticking probability becomes unimportant onc e evaporation be 



comes dominant. We assume the same set of adsorption energies (-E a ds) as lGarrod fc Herbst 
( 120061 ) . corresponding to a grain surface covered by water ice. The adsorption energies 



for deuterated species are set to the values used for their normal isotopes, with the ex- 
ception of D atoms, whose adsorption energy is set 21 K higher (471 K) than that of 
H atoms (450 K), following ICaselli et al.l ( 120021 ). The grain-surface species can desorb to 
the g as phase via thermal e vaporation and two non-thermal processes: cosmic-ray desorp- 
tion (IHasegawa et al.l 119931) and desorption via exothermic associati on reactions on grain 



surfaces ( jWilliams 



1968 



Garrod fc Herbstl l2006t iGarrod et al.l 120071 ). Recent laboratory 



work an d molecular dynamics simulations i ndicate that photo de sorption could also be im 



Arasa et al. 



portant (lAndersson fc van Dishoeck 1 120081 ; lOberg et al.l l2009al Jbl; iMunoz Caro et al. 1 12010 



20111 ) . There are several possible mechanisms for photodesorption: Molecules 



are dissociated into radicals and atoms, which then desorb separately to the gas phase; the 

disso ciation products may first recombine on the grain and then desorb ( lAndersson &: van Dishoeck 
20081 ): or the dissociation products may kick out neighboring ice-mantle species. Our model 
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includes the second of these mechanisms, through the explicit inclusion of two separate 
processes; the photodissociation into radicals of grain-surface species by interstellar and cos- 
mic ray- induced UV photons (see below), and t he partial evaporation of the products of 
exothermic surface reactions ( IGarrod et al.l 120071 ). 



We assume the Langmuir-Hinshelwood mechanism for grain surface reactions; species 
can diffuse on grain surfaces by thermal hopping and react with each other when they 
meet. No quantum tunneling was assumed in the migration, even for H and D atoms 
( IGarrod et al.l 120071 ) . The barrier for migration is set at 50% of the adsorption energy. We 
adopt the modified-rat e approach for grai n-surface reactions that include H and D atoms 
( IStantcheva et al.ll200ll ; ICaselli et al.l 120021 ) . In paper I, we did not adopt the modified-rate 
approach, since its significance is reduced in the model when assuming only thermal hopping 
of H atoms on the gra in surface, compared with a model that includes quantum tunneling 
( Muffle fc Herbstl 120001 ) . In the present work, however, the grain-surface reactions with D 
atoms could still be accretion-limited rather than migration-limited, and thus the modified 
rates can be important. 

Extinction of interstellar UV radiation is calculated from the column density of hydrogen 
nuclei (A%) from the core outer edge to each shell by the formulation Ay = ./Vh/(1.59 x 
10 21 cm 2 ) mag. Initially, the visual extinction from the outer edge (r = 4x 10 4 AU) to the core 
center is about 5.5 mag, and the outermost shells in which we calculate chemistry are located 
at 1 x 10 4 AU, where A v ~ 1 mag. Assuming that our model core is embedded in ambient 
clouds, we add 3 mag to the visual extinction obtained above, and ignore photodissociation of 
CO and H 2 , which should be self- and mutually-shielded from the interstellar UV at A v > 1 
mag (e.g. Ivan Dishoeck fc Blacklll988l ). We also neglect the photodissociation of HD and D 2 , 
since it is not clear if and how much they are shielded by H 2 . This assumption would not 
much affect our results for molecular D/H ratios, since the interstellar UV is extinguished 
by > 3 mag. The rates of photodissociation and photoionization by cosmic ray-induced 
UV radiation, on the other hand, are roughly independent of A v and proportional to the 
cosmic-ray ionization rate, although the coefficients vary among molecular species. The 
photodissociation rates of HD and D 2 by cosmic-ray induced UV are set to be the same as 
for H 2 . UV radiation from the central protostar is not considered in our model; since the 
model is ID, UV radiation is absorbed in close vicinity to the source. 
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3. Results 

3.1. Temporal Variation of Molecular Distributions 

Figured {a — b) shows the radial distribution of simple molecules at t = —5.6 x 10 2 yr, 
4.3 x 10 2 yr, and 9.3 x 10 4 yr. These simple molecules comprise four major neutrals: CO, 
NH3, H2O and N2. Major ions and electron abundance (i.e., the fractional ionization) are 
also plotted. A simple view of the radial distribution of major species is as follows: as the 
radius gets smaller, and the density higher, the species first tend to be depleted from the 
gas onto the ice. Finally, as the temperature rises, sublimation occurs and the species return 
rapidly and sharply to the gas at the so-called sublimation radius. Little change is seen in 
gas-phase values at yet smaller radii, due to the rapid collapse. Such a radial distribution of 
gas-phase abundance is consistent with the "drop ab undance profile" derived from multi-lin e 



observations of CO and H 2 toward protostars (e.g. iJeirgensen 112004 : ICoutens et al. II2012I ). 
The sublimation radii of CO, NH3, H2O and N2 increase with time, as the core temperature 
rises, so that the gas-phase material extends farther from the core. The CO depletion factor 
is defined as the ratio of the canonical value of CO, which is close to the gaseous elemental 
abundance of carbon (8 x 10 -5 in our model), to the actual abundance of gaseous CO. The 
depletion factor increases inward as the density increases, until it suddenly decreases at the 
sublimation radius. The peak value of the CO depletion factor decreases with time, because 
the density decreases. We note that the CO depletion factor does not reach unity even 
inside the sublimation radius in our model; i.e. CO abundance inside the sublimation radius 
is smaller than the canonical value, especially at early stages, because CH 3 OH and CH 4 are 
slightly more abundant than CO. Low CO abunda nces inside the subl i mation radius have 



recently been observed t oward several protostars (jAlonso-Albi et al. I l2010l : lYildiz et al. 



20121 : iFuente et al. 1120121 ) 



The fractional ionization decreases toward the center, because it is generally proportional 
to n H . It is noteworthy that HCO + , which is often used as a radio tracer of dense gas, 
decreases as two species with a higher proton affinity than CO (H2CO and NH3) sublime 
into the gas phase. 

The distribution of complex organic species is shown in Figured (c). At t = —5.6 x 10 2 
yr, the temperature is lower than ~ 10 K at > 100 AU (Figured]). In these cold regions, 
CH3OH ice is as abundant as 1.62 x 10~ 5 , and HCOOH ice and CH 3 CN ice are ~ 10~ 9 
relative to hydrogen nuclei. CH3OH ice is formed by the hydrogenation of CO ice, while 
HCOOH is mainly formed in the gas phase via the dissociative recombination of C^O^, 
which is a product of the radiative-association reaction HCO + + H 2 0. A fraction of the 
gas-phase HCOOH is adsorbed onto grains before subsequent gas-phase reactions. CH 3 CN 
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is formed by association reactions (which add C and H to CN) on grain surfaces. Because 
these species are formed by low-temperature chemistry, the abundances at this stage depend 
on the duration of the pre-collapse model assumed. HCOOCH 3 ice, on the other hand, 
increases inwards at 10 AU; it is formed at warm temperatures, starting at 30-40 K. 

At t = 4.3 x 10 2 yr, the warm region of T > 100 K expands to ~ 10 AU, where complex 
organic species sublime. At t = 9.3 x 10 4 yr, many organic species are forme d via grain surface 



reacti ons and gas phase reactions of sublimated species at T > 30 K (IGarrod fc Herbst 



( 120061 ); Paper I). Our model shows that the abundance of large organic species increases with 
time at the central region not only because the dust temperature exceeds their sublimation 
temperatures, but also because the duration of the warm chemistry is longer for fluid parcels 
in later stages than those in earlier stages (see Figured]). It should also be noted that the 
abundance rise is not homogeneous; relative abundances among complex organic species (e.g. 
CH3CHO/HCOOCH3) vary with time. 

Figure EJ^cT) shows the distribution of carbon-chain species. At t — 9.3 x 10 4 yr, carbon 
chains increase inwards at the sublimation radius of CH 4 (T ~ 25 K, r ~ 2000 AU), showing 
that the chemistry is indeed described by WCCC; the CH 4 reacts with C + to form C2H3 , 
and unsaturated hydrocarbons are formed yia subsequent gas -phase re actions and a parti al 
contribution of grain surface hydrogenation (ISakai et al.l (120081 ) ;Paper I; IHassel et al.l ( 120081 )). 
At earlier stages (t = —5.6 x 10 2 yr and 4.3 x 10 2 yr), the formation of carbon-chains at the 
CH 4 sublimation radius is much less significant. Although the CH 4 abundance is as high 
as 2.5 x 10~ 5 , the other reactant to start WCCC, C + , is deficient. Specifically, while the 
abundance of C + at the sublimation radius of CH 4 is as high as ~ 1 x 10~ 10 at t = 9.3 x 10 4 
yr, it is 4 x 10 -15 at t = —5.6 x 10 2 yr, because of the high density. Although some carbon 
chains are formed at the CH 4 sublimation zone at t = 4.3 x 10 2 yr, the increment is much 
smaller than that at t = 9.3 x 10 4 yr. The abundances of carbon chains (e.g. C3H2) at the 
central region of these early phases are thus determined mostly by the sublimation of icy 
carbon chains formed in the cold era. 

It should also be noted that C 4 H and C2H are depleted at radii less t han a few hundre d 
AU at t = 9.3 x 10 4 yr, which is consistent with observations of L1527 ( ISakai et al.ll2010l ). 
In these hot- temperature region s they react with H 2 and are converted to C 4 H 2 and C2H2, 



respectively (IHassel et a 



region, while ISakai et al 



2011 ). C3H2, on the other hand, does not deplete at the central 



(2010) found that the cyclic form, C-C3H2, does indeed exhibit a 



central dip. IHassel et al.l ( 1201 ll ). who distinguish C-C3H2 from H 2 CCC, also show that c- 
C3H2 survives in hot regions (~ 100 — 200 K) for several 10 5 yr, while the carbene form is 
destroyed via reaction with H2 over a shorter time scale at T > 140 K. Additional destruction 
paths for C-C3H2 are therefore required to account for the observations. 
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3.2. Molecular D/H ratios 

Figure [3] shows the spatial distribution of selected mono-deuterated species, as well as 
D^ , and their normal isotope counterparts at t — —5.6 x 10 2 yr, 4.3 x 10 2 yr and 9.3 x 10 4 yr. 
The D /H ratios of assorted molecules are shown in the bottom panels in Figure |3j In low- 
temperature chemistry (T ~ 10 K), gas-phase molecules are enriched with deuterium due to 
exothermic exchange reactions such as H^ + HD — > H 2 D + + H2 (and analogous reactions 
for multiply deuterated species; Millar, Bennett & Herbst 1989; Roberts, Herbst & Millar 
2004). The backward reactions, e.g. H2D + + H2 — > H3 + HD, are endothermic (typically 
by several 100 K) and thus are negligible at low temperature (~ 10 K). The D/H ratio for 
H3 is not governed by these two reactions alone, as H3 and its deuterated isotopomers are 
destroyed by reactions with CO and electrons. In the zone where gaseous CO is depleted, 
however, the D/H ratio for H^ is further enhanced. For example, at t — —5.6 x 10 2 yr, 
Djj~ becomes the major cation in the CO depletion zone at radii of a few tens to several 
hundred AU (see Fig. [3D- Although the exothermic exchange reactions are limited to several 
species, the high D/H ratio propagates to other molecules via gas-phase reactions; e.g. H 2 D + 
gives its deuteron to other molecules and atoms via ion-molecule reactions. It is essential 
to include multi-deuterated species, such as D3 , to calculate the molecular D/H ratio, since 
the multiply-deuterated species can propagate D to other species more efficiently than can 
mono-deuterated ones. The high D/H ratio also propagates to the grain surface; deuterated 
H^ dissociatively recombines with electrons to produce D atoms, which are adsorbed onto 
grains and deuterate the ice-mantle species. 

Once the core gets warmer than ~ 20 K, CO sublimation terminates the extreme deuter- 
ation (XD/XH> 1) of H^ . At the same time, the endothermic exchange reaction, H 2 D + + 
H 2 — > H3 + HD, becomes effective and dominates in the destruction of deuterated H^ at 
higher temperatures. Since H 2 is very abundant, the ratio of H 2 D + to H^ rapidly reaches 
its equilibrium value, which is ~ 1.3 x l Q- 5 fr/300K)- a8 exp(-230K/T ), i.e. 4 x 10" 3 and 



3 x 10" 4 at 30 K and 50 K, respectively (Millar. Bennet & Herbstlll989f ). Hence H 2 D+ and 



HD^ may be good observational temperature probes in low- to intermediate-temperature 
regions (< 50 K) of the core. 

Neutral species have high D/H ratios even in warm regions. These neutral species, such 
as NH 3 , CH 4 , and CH 3 OH, are formed more efficiently via grain-surface reactions than by 
gas-phase reactions in the cold era. Inside the sublimation radii, the abundances of the 
sublimates overwhelm those of the products of gas-phase reactions. The icy material is 
highly deuterated by surface reactions with D atoms as described above, although the D/H 
ratios inside the sublimation radii are lower than the ratios of their gaseous counterparts 
in the CO depletion zone, since ices are formed and accumulated throughout the cold era, 



13 



including the time before significant CO depletion occurs. For the same reason, the D/H 
ratios of the sublimates (i.e. inside the sublimation radii) are higher in earlier stages only 
by a factor of a few. The D/H ratio of these gaseous neutral species is constant down to 
r = 2.5 AU, since their destruction time scales are longer than a free-fall time scale. 

At t = 9.3 x 10 4 yr, various complex organics are formed via both gas-phase and grain- 
surface reactions; various carbon chains and large organics are formed via WCCC at T > 25 
K and via radical-radical surface reactions at T > 30 K, respectively. Their D/H ratios are 
high, mainly because they are formed from molecules with high D/H ratios. As an example, 
methyl formate (HCOOCH3) can be formed by the surface radical-radical reaction between 
the radicals HCO and CH 3 0. Since both HCO and CH 3 have high D/H ratios, DCOOCH3 
and HCOOCH 2 D are also formed. The HCOOCH 2 D is three-times more abundant than 
DCOOCH3, since we assume the statistical branching rather than following the position of 
D in each reaction. 

It is well known that the exchange reactions of CH^ + HD and C^H^ + HD are more 
exothermic than H^ + HD, and thus can be effective even at relatively warm temperatures 
(several tens of K). However, they do not much contribute to the deuteration of WCCC 
species in our model, probably because carbon chains and CH4 already have high D/H 
ratios when the fluid parcels enter the warm region, and because the duration of the warm 
temperature (around 30 K) is short. The D/H ratios of carbon-chain species change little, 
even if we neglect these two exchange reactions. 

The exchange reaction of OH + D — > OD + H, on the other hand, has an even higher 
exothermicity (AE/k ~ 810 K), and significantly deuterates OH at t — 9.3 x 10 4 yr at radii 
of 100 - 1000 AU, where HCOOH has a higher D/H ratio than other organic species due to 
the reaction of sublimated H 2 CO with OD. 

Radial distributions of assorted multiply-deuterated species are shown in Figure HI to- 
gether with mono-deuterated and normal isotopologues. The abundance ratios of multi- to 
mono-deuterated species are mainly determined statistically, e.g. XD 2 /XHD ~ XHD/XH 2 . 
For CH3OH and H 2 CO, however, the multi- to mono ratio is lower than the mono-deuterated 
to normal isotope ratio. This effect could be due to the activation barriers in their forma- 
tion path (i.e. in the reactions H + CO and H + H 2 CO); the tunneling rate through the 
activation barrier is lower for heavier isotopomers (see discussion in §4.3). For NH3 and 
CH 4 , on the other hand, the multi-deuterated species are more abundant than expected 
from their mono-deuterated to normal isotope ratio. Although they are mainly formed on 
grain surfaces, gas-phase formation also contributes. In the gas-phase, multiply-deuterated 
NH3 and CH4 are efficiently formed in the CO depletion zone. Since the abundances of 
multiply-deuterated species are generally low, the efficient formation in the CO depletion 
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zone effectively increases its abundance relative to the normal isotope. At t — 9.3 x 10 4 yr, 
multi-deuterated NH 3 and CH 4 are much less abundant than in early stages, because the 
CO depletion is much less significant. 

In our model, we do not consider the layered structure of ice mantles. In reality, ice 
mantles have layered structures and the ices with the highest D/H ratios, which are formed 
at the latest stage of the cold prestellar phase, would lie on the very surface layer of the ice 
mantle, and would be the first to be sublimated. This effect creates radial gradients of D/H 
ratios. For example, the CH 4 formed in the CO depletion zone would have a higher D/H 
ratio than that formed in the early molecular cloud stage. Since the CH4 formed in the early 
stage would coexist within H 2 ice, mainly the highly deuterated CH 4 (in the surface layers) 
would be sublimated at T ~ 25 K; e.g. at t — 4.3 x 10 2 yr the gaseous CH 4 abundance would 
be lower and CH3D/CH4 ratio would higher in r = 10 — 100 AU than shown in Figure [3J 
Once the temperature reaches T > 150 K, however, the dominant ice component, H2O ice, 
sublimates together with other species formed in the early stages of molecular clouds. Then 
the abundance and D/H ratio of CH4, for example, would be the same as shown in Figure [3] 
inside ~ 10 AU at t = 4.3 x 10 2 yr. 



4. Discussion 

4.1. Effect of 2D structure 

So far, we have adopted a spherically-symmetric model of star-forming cores. In the 
spherical model, fluid parcels fall to the central star over a short time scale. But in reality, 
dense cores are rotating with uj ~ 10~ 14 s _1 , and the spherical symmetry is broken in the 
central regions as the collapse proceeds. Because of angular-momentum conservation, the 
centrifugal force increases for infalling fluid parcels and eventually balances the gravity at 
the centrifugal radius, which is given by the expression 



GM V°- 07 P C / V lxl 0~ 14 s-V \M 

to form a circumstellar disk, where r and M are the size and mass of the initial dense core, 
respectively. The result is that the fluid parcels stay in warm dense regions for a longer time 
scale than occurs in the spherical model. 

In order to investigate the chemistry in this prolonged warm phase, we have performed 
a pseudo-time-dependent model (i.e., fixed density and temperature). For the initial abun- 
dances, we adopt the molecular abundances at 30 AU in the final step of the spherical 
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collapse model (t = 9.3 x 10 4 yr), considering that the typical centrifugal radius is ~ 100 AU 
and that molecular abundances at 2.5 - 100 AU are almost constant in the collapse model 
at 9.3 x 10 4 yr. All ices have already evaporated to the gas phase in this initial abundances 
(Figure [2]). 

Figure [5] shows the temporal variation of molecular abundances at T = 40 K (a — h) and 
T = 150 K (c — d). We choose these temperatures because grain-surface reactions would be 
active at T = 40 K, while gas-phase reactions would dominate at T = 150 K. The density is 
set to be nn = 1.487 x 10 8 cm -3 , which is the density at r = 30 AU in the ID model. These 
parameter choices are rather arbitrary, but there should be regions with these densities and 
temperatures in the disks, because a real disk has density and temperature gradients in both 
the vertical and radial directions. We simply assume that the infalling fluid parcel reaches 
and stays in such regions in the disk. We also performed models with higher (xlO) and lower 
(xO.l) densities, but the results are qualitatively the same; the timescale of the che mical 
evolution described below depends only slightly on the density (INomura et al. 1 120091 ) . At 
T = 40 K, C2H6 increases and becomes the dominant carbon chain species. It is formed via 
grain-surface hydrogenation of C 2 H 4 , which is produced by the g as-phase reaction of CH + 
CH 4 . It should be noted that C 2 H 4 (£ ads = 3487 K in our model; iGarrod fc Herbstl (bood )) 
can marginally freeze onto grains, while CH (E^s = 925 K) and CH4 (E a a s = 1300 K) are 
easily desorbed at 40 K; adsorption of C2H4 onto the grain surface works as a sink in the 
gas-phase chemical network. Similarly, C7H4 ice (-E a ds — 7487 K) is the dominant carbon 
chain at t = 10 6 yr at 150 K. Since the sublimation temperature of C7H4 is close to 150 
K, it depletes onto grains before additional carbon-chain growth in the gas phase at high 
densities, while hydrogenation on grain surface is inefficient at this high temperature. 

The D/H ratios of the carbon chains remain high at T = 40 K, since the mother 
molecules (e.g. CH 4 ) are highly deuterated in the cold phase. Also, the exothermic exchange 
reactions CHj~ + HD and C 2 H^ + HD, which have higher exothermicities than + HD, 
contribute to keep the D/H ratios high. At T = 150 K, on the other hand, the D/H ratio of 
carbon-chain species decreases gradually, because the endothermic exchange reactions (e.g. 
CH 2 D + + H 2 and C 2 HD + + H 2 ) can now proceed efficiently It is interesting that species that 
undergo direct exchange reactions with HD, such as H3" and C 2 H 2 f , work as a source of D/H 
enhancement in the network at low temperatures, but easily lose their own D enhancements 
at high temperature. 

Among complex organic species, the HCOOCH3 abundance increases from 1.6 x 10 -9 
(in the gas phase) to ~ 10~ 7 (in the ice) via the grain-surface reaction of HCO + CH 3 in a 
few times 10 yr at T = 40 K, while the abundances of other complex species do not change 
significantly within this timescale. At T = 150 K, on the other hand, a significant amount 



- 16 - 



of CH3OCH3 is formed in the gas phase from CH3OH in ~ 10 5 yr. In spite of the formation 
at high temperatures, the D/H ratios of CH3OCH3 is high, because the mother molecule, 
CH3OH, is highly deuterated in the cold phase. At t > 10 5 yr, complex organic species are 
destroyed via protonation and subsequent dissociative recombinations. 

Another important feature missing in our spherical model is outflows. It creates a 
cavity in the envelope, through which the UV and X-ray from the protostar may escape 



Photodissociation and photoionization play important roles near the cavity wall ( jVisser et al. 



201ll ). or the surface of the circumstellar disks, if the envelope mass is small. The fraction 
of the mass going through such PDR-like regions is, however, rather limited: < several % 



(jVisser et al.ll201ll ). Our ID model would be a reasonable approximation for the rest of the 



envelope regions, relocated from the outflow cavity. 



4.2. Comparison with observations 

Table [1] lists estimated abundances with respect to H2 from observations of oxygen- 
containing organic species in the protostellar core IRAS 16293-2422. We distinguish the 
observations by single-dish telescopes and interferometers. IRAS 16293-2422 is a binary sys- 
tem; interferometric observations distinguish the two sources IRAS16293A and IRAS16293B, 
while single-dish observations integrate the emission from the whole system. As can be seen 
from the diverse results, it is not straightforward to estimate molecular abundances in pro- 
tostellar cores, because the cores have temperature and density gradients along the line of 
sight, and because the high-temperature regions, in which these complex organic species sub- 
lime, are mostly spatially unresolved even with the interferometer observations. Abundance 
estimates therefore vary significantly depending on the assumed temperature and density 
structure of the core. 

It is noteworthy that interferometric observations tend to derive lower ab undances of 



l arge c omplex molecules than single dish observations corrected for beam dilution. lOberg et al 



(120111 ) also noted such trends, and pointed out that cold extended component contributes 
significantly to the single dish observation. Thanks to the higher spatial resolution, inter- 
ferometric observations could derive more reliable abundances in the central hot (> 100 K) 
regions. 

Table [T] also lists the abundances in our model near the center (r = 30 AU) of the 
protostellar core at t = 9.3 x 10 4 yr and in our pseudo-time-dependent model with T = 40 
and 150 K at t = 10 5 yr from the start of the "disk". The core radius of 30 AU is arbitrary, 
because the molecular abundances inside ~ 100 AU are almost constant due to rapid infall; 
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it simply represents the abundances at T > 100 K. As for the pseuo-time-dependent model 
at T = 40 K, we list the ice abundances of complex organic species, assuming that they 
eventually accrete to smaller radii to be sublimated to the gas phase. As can be seen, the 
abundances of H 2 CO and CH 3 OH are higher in our core model than derived from obser- 
vation. In the pseudo-time-dependent models, on the other hand, the H 2 CO abundance is 
far too low at both temperatures at the time used, while the methanol abundance for the 
150 K model is in good agreement with observation. Note though that the abundances of 
H 2 CO and CH 3 OH are time-dependent; these species decrease with time at 150 K (Figure 
[5]). We note that the spatial e xtent of the CH3OH sublimation zone (i.e. T > 100 K) in the 



protostellar core is ~ 100 AU flSchoier et aubooi l. which is comparable to the typical cen- 



trifugal radius of dense cores. Therefore CH3OH can be partially destroyed in the prolonged 
warm temperature caused by rotational support. It is noteworthy that a recent ALMA ob- 
servation of methyl formate revealed a velocity gradient which is consistent with rotating 
gas around IRAS 16293 A, although it is not ruled out that the velocity gradi ent could orig- 



inate in two unresolved velocity components such as occur in a close binary (IPineda et al. 



20121 ). Another explanation for the disagreement of CH3OH and H2CO abundances between 
our core model and the observation might be that these abundances are underestimated by 
observations. Observational constraints of molecule abundances in the warmer regions are 
uncertain because of beam dilution effects, the presence of multiple systems, and the possible 
presence of cavities and outflows. Indeed, interferometric observations which distinguish the 
two protostars in IRAS16293 system gives higher abundances of CH3OH than single dish 
observations. 

Concerning the other molecules in the table, the model abundances for HCOOCH3, 
HCOOH, and CH3CN in the core model are within the range of the estimated abundance 
from the observations. CH3OCH3, on the other hand, is significantly underestimated in 
our core model. In the pseudo-time-dependent model at 150 K, however, the CH3OCH3 
abundance increases with time and agrees with the observation at t ~ 10 5 yr. It is formed 
by the recombination of CH3OCH4 in the gas phase. We note that the branching ratio 
for producing CH3OCH3 in the recombination is only 5 % in our model. If we assume a 
higher branching ratio, CH3OCH3 becomes more abundant. Abundances of other species 
in the pseudo-time dependent model at 10 5 yr are also consistent with observations, except 
for H2CO. This may indicate that the gases in the central hot region of IRAS16293 may 
have experienced prolonged warm temperature chemistry compared with the spherical infall 
model due to partial rotational support. Observations with higher angular resolution are 
highly desirable for more detailed comparison. 

Our co r e and pseudo-time-dependent models can also be applied to carbon-chain species. 



Sakai et all (120 101 ) observed carbon chains towards LI 52 7, and showed that the intensity 
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distribution of C3H2 (43,2 — 42,3) is well reproduced, if its abundance jumps from 2.7 x 10~ 9 to 
2.7 x 10~ 8 in the temperature range T = 23 — 29 K. These peak abundances and temperature 
ranges are in reaso nable agreement with our core model at t — 9.3 x 10 4 yr and r ~ 10 3 AU. 



Sakai et all (120101 ) found that C 2 H, C 4 H, and C 3 H 2 have a slight dip at the central region. 
In our spherical model, the central holes of the C2H and C4H abundances are reproduced, 
but C3H2 increases inwards, as discussed earlier. In the pseudo-time-dependent model at 40 
K, C3H2 is quickly depleted onto the grains, just like gaseous C 2 H 6 in Figure [51 The central 
dip of C3H2 thus might be reproduced, if the central region is dominated by the cold disk 
material rather than a warm infalling envelope. It is possible that the disk midplane is much 
cooler than the envelope because of the higher optical depth. 

Let us now move on to deuterated species. Table [2] lists the estimated molecular D/H 
ratios in various protostellar cores and our model D /H ratios. As in Table [Tj we list model 
D/H ratios at t = 9.3 x 10 4 yr of the protostellar core at r = 30 AU, except for C4H, and 
in our pseudo-time-dependent models. For the C4D/C4H ratio in the core model, we list 
the ratio at r = 250 AU, where C4H abundance reaches its maximum. Our D/H ratios for 
the carbon chain species C4H and C3H2, and for NH 3 , are in reasonable agreement with the 
observations of L1527. On the other hand, the observed D/H ratios of NH 3 in NGC1333 
and H2CO in several protostellar cores are much higher than in our models (either core or 
disk). It should be noted, however, that in our core model these ratios are higher in the 
central warm region at early stages (e.g. t = 4.3 x 10 2 yr) and in the CO depletion zone 
than listed in Table 3 (Figure HJ). The model D/H ratios vary significantly depending on 
the evolutionary stage of the core, and the size of the warm regions, while the observational 
D/H ratios thus depend on beam size (whether the beam traces mainly the cold component 
or the central warm regions). Such dependence also explains the variety of observed D/H 
ratios towards protostars. 

The CH 2 DOH/CH 3 OH ratio reaches a high of 0.65 in the observations, while it is only a 
few per cent in our models. The comparison is even worse for CHD2OH/CH3OH. Unlike NH3 
and H 2 CO, CH 3 OH is not formed in the gas phase, and its gaseous abundance is very low in 
the CO depletion zone. Hence it is unlikely that the high CH 2 DOH/CH 3 OH ratio originates 
in the CO depletion zone. In our model, CH 2 DOH is formed by the addition of D during 
the hydrogenation of CO on grain surfaces; we have not included substitution reactions of 
H2CO + D , or the abstraction reaction CH3OH + D , whic h are observed in laboratory ex- 



perim ents (INagaoka et al.l 120051 120071 ; iHidaka et al. 1 120091 ) . iTaquet. Ceccarelli fc Kahane 



(120121 ) recently performed gas-grain chemical models at constant density and temperatures, 
and showed that these reactions are crucial to reproduce the observed D/H ratios of H 2 CO 
and CH3OH. In the following subsection (§4.3) , we add the substitution and abstraction 
reactions on grain surfaces to our network to see how strongly the CH2DOH/CH3OH abun- 
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dance ratio is enhanced. 



The D/H ratio of water ice in the envelope at t — 9.3 x 10 4 yr is ~ 1 %, which is 



marg inally consistent with the upper limits obtained towards low-mass protostars (IParise et al 
20031 ). Recently gaseous HDO has also been intensively searched for towards protostellar 
cores using ground-based single-dish telescopes and the Herschel Space Observatory. Ta- 
ble |3] lists the estimated HDO abundance and HDO/H 2 ratio in the gas phase in the hot 
corino regio n (T > 100 K) and outer envelope (T < 100 K) towards NGC 1333-IRAS2A and 
IRAS16293 flLiu et al. IlioiH ICoutens et al. Il2012h ~ Towards IRAS16293, an absorbing layer 
is added in front of the envelope to account for the observed absorption lines. Our model 
HDO/H2O ratio in the hot region is consistent with the observations, while our HDO abun- 
dance is higher by an order of magnitude (see panel (a) of t = 9.3 x 10 4 yr in Figure 3). We 
note that the HDO (and H 2 0) abundance could be underestimated in observations, just like 
CH3OH and H 2 CO, because of the complex (and poorly-constrained) structures of the source 
such as the presence of multiple systems (IRAS 16293), and the possible presence of cavities. 



outflo ws and disks. In addition to the targets listed in Table |3j IjQrgensen fc van Dishoeck 



(120101 ) searched for gaseous HDO towards NGC1333-IRAS4B using the Sub-Millimeter Ar- 
ray (SMA). The line was not detected, which constrains the HDO/H 2 ratio < 6 x 10~ 4 
in the central hot region (~ 50 AU). This upper limit is much lower than predicted in our 
model. Since gaseous HDO in the central region originate in ice formed in the cold era, 
the variation of the HDO/H 2 ratio among protostars could originates in t he temperature 



varia tion in the very early stages of molecular clouds or translucent clouds ( ICazaux et al. 
201 lh . 



In the outer envelope, our HDO abundance and HDO/H 2 ratio at > 10 3 AU at t — 
9.3 x 10 4 yr are in reasonable agreement with NGC1333-IRAS2A and absorbing layer of 
IRAS16293 (Figure 3). But the observed values in the envelope of IRAS16293 is lower than 
in our model at t = 9.3 x 10 4 yr. Since the molecular column densities in the envelope 
are lower in earlier stages, a model at slightly earlier stage than t = 9.3 x 10 4 yr might be 
more consistent with the envelope HDO abundance in IRAS 16293. It should also be noted 
that in our model at t = 9.3 x 10 4 yr, HDO/H 2 ratio is close to unity around 500 AU, 
which is a part of the outer envelope (Figure E]). This high HDO/H 2 ratio is caused by the 
exothermic exchange reaction o f OH + D; the rate coefficient is set to k = 1.3 x 10~ 10 cm 3 s~ 1 
(IMillar. Bennet &: Herbstl 119891 ). which is consistent with a detailed theoretical calculation 
bv lAtahan et al. I J2005F fl.5 x 10" 10 cnrV 1 at 50 K). Such local HDO/H 2 enhancement 
has not been observed thus far. Although the effect of the OH + D reaction seems to be 
transient and local in our model, a further search for such enhancements could constrain 
O-chemistry in the protostellar envelope. 
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Table 1: Observed and Model Gas-Phase Molecular Abundances in IRAS 16293-2422 



Species 




IRAS 16293-2422 




model a 






single dish 


interferometer 


core 


40 K 


150 K 


H 2 CO 


1.0(-7) b 


l.l(-7) c 


l.l(-5) 


5.4(-10) d 


l.l(-ll) 


CH 3 OH 


1.0(-7) c 


9.4(-8) c , 6.8(-7) f , 3.1(-7) s 


3.6(-6) 


2.9(-6) d 


5.2(-7) 


HCOOCH3 


4.0(-7) h 


7.5 (-9) i ,4.3(-9) f , 2.6(-9) s 


3.1(-9) 


1.7(-7) d 


8.8(-10) 






4.9 (-IO)J ,> 1.2(-8) k 








HCOOH 


6.2(-8) h 


2.5(-9) k 


l.l(-8) 


3.3(-8) d 


1.2(-8) 


CH3OCH3 


2.4(-7) h 


7.6(-8) c 


1.2(-10) 


3.9(-10) d 


3.3(-8) 


CH 3 CN 


7.5(-9) ! , 1.0( 


-8) h 4.8 (-lO) 1 , 2.3(-ll) j 


5.6(-9) 


1.3(-8) d 


4.3(-9) 



"Molecular abundances at r — 30 AU and t — 9.3 x 10 4 yr in the protostellar core model, and in the 
pseudo-time-dependent model of T = 40 K and 150 K at t = 1 x 10 5 from the start of the "disk". 



1 Maret et al. 


(2004) 


c Chandler et al. ( 


2005) 



d ice abundance 
Maret et all |2005[) 



^ iKuan et al. 1 (|2004l ), IRAS16293A 
4Kuan et al.1 |2004l), IRAS16293B 
^Cazaux et all |2003l ) 



4 IRAS16293 A (|Bottinelli et al.ll2004l) 
J 'IRAS16293 B |Bottinelli et al.ll2004l ) 
feemiiian fc Hollis I |2006l ) 



Schoier et al.l (2002) 



Table 2: Molecular D/H ratios in protostellar cores 



Species 


D/H ratio 


source 


core 


model 3, 
40 K 


150 K 


reference 




NH 2 D/NH 3 


0.28 


NGC1333 


0.015 


0.014 b 


0.012 


Roueff et al. 


(2005) 




0.04 


L1527 








Sakai et al. 1 


'2009) 


ND3/NH3 


2.1 x 10~ 4 


NGC1333 


1.6 x 10~ 6 


1.5 x 10~ 6b 


5.1 x 10~ 7 


Roueff et al. 


(2005) 


HDCO/H 2 CO 


0.094-1.7 


7 protostars 


0.026 


0.013 b 


_c 


Parise et al. 


(2006) 


D 2 CO/H 2 CO 


0.046-0.44 


7 protostars 


3.3 x 10~ 5 


2.9 x 10~ 5b 


_c 


Parise et al. 


(2006) 


CH2DOH/CH3OH 


0.37-0.65 


4 protostars 


0.024 


0.026 b 


0.023 


Parise et al. 


(2006) 




< 0.030 


L1527 








Sakai et al. 1 


;2009) 


CHD2OH/CH3OH 


0.074-0.25 


4 protostars 


1.7 x 10~ 4 


2.0 x 10~ 4b 


1.5 x 10~ 4 


Parise et al. 


(2006) 


C4D/C4H 


0.018 


L1527 


0.033 d 


c 


_c 


Sakai et al. ( 


'2009) 


C 3 HD/C 3 H 2 


0.071 


L1527 


0.033 


0.016 b 


0.0083 


Sakai et al. 1 


;2009) 



"Molecular abundances at r — 30 AU and t — 9.3 x 10 4 yr in the protostellar core model, and in the 
pseudo-time-dependent model of T — 40 K and 150 K at t = 1 x 10 5 from the start of the "disk". 

b ice abundance 

C D/H ratio is not listed because the abundance of the molecule is too small compared with observations. 
d D/H ratio at r = 250 AU, where C4H abundance reaches its maximum. 



Table 3: HDO abundance and HDO/H20 ratio in the gas phase in protostellar cores 



regions reference 
hot (T > 100 K) envelope (T < 100 K) absorbing 



NGC1333-IRAS2A Liu et al. (2011) 

HDO abundance 6.6 x 10~ 8 - 1.0 x 10" 7 9 x 10~ n - 1.8 x 10~ 9 

HDO/H2O > 1 % 0.9 - 18 % 

IRAS16293-2422 Coutens et al. (2012) 

HDO abundance 1.7 x 10~ 7 8 x 10" 11 

HDO/H 2 1.4 - 5.8 % 0.2 - 2.2 % 4.8 % 
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4.3. Deuteration of CH 3 OH on grain surfaces 



The grain-surface react ion network for H 2 CO and C H3 OH derived from labo ratory ex- 
periments is summarized in IWatanabe &: Kouchi I (120081 ) and iHidaka et al. I (j2009[ ) . Most of 
the reactions are already included in our original gas-grain model, except for the substitution 
H 2 CO + D HDCO + H, abstraction CH 3 OH + D CH 2 OH + HD, and their deuterated 
analogues (Figure EJ). Here we add these reactions to our model. 

Laboratory experiments indicate that the abstraction reaction CH3OH + D produces 
CH 2 OH but not CH 3 0, so the subsequent hydrogenation reaction can form CH 2 DOH but not 
CH3OD, which is i n agreement with the observed low abun dance of CH„D 3 _ n OD compared 
with CH„D 3 _ n OH ( Parise et al. 112004 iNagaoka et al.ll2005l ). In our original deuterated net- 
work model, however, we did not follow the position of the deuterium atoms in the reactions, 
and we did not fully distinguish the various possible isomers of radicals, such as CH 2 DO, 
CHDOH and CH 2 OD, which are produced in the sequence of deuterated methanol forma- 
tion. Our model thus underestimates (dilutes) the isomer ratio CH n D 3 _ n OH/CH n D 2 _ n OD, 
although we distinguish isomers in the abstraction reactions if they are in our original species 
list. In spite of the imperfect distinction of the isomers, CH 2 DOH is more abundant than 
CH3OD in our model. 

The activation barriers of the substitution and abstraction reactions have recently been 
investigated using theoretical calculations (e.g. iGoumans 1 1201 lh : however, the rate is deter- 



mined not only by activation barrier height. Since the reactions proceed via tunneling at 
low temperatures, the tunneling mass and the shape (width) of the potential barrier are also 
important. Thus, the accurate evaluation of reaction rates from theoretical calculations is 
challenging. Rather than incorporate activation barriers from theoretical calculations, here 
we utilize the rates derive d from the laboratory experiments by IWatanabe fc Kouchi I (120081 ) 
and lHidaka et al. I (120091 ) . who estimated relative rate coefficients of various reactions to the 
rate of H + CO. In our model, rate coefficients for grain-surface reactions are proportional 
to a frequency factor times exp [—2w/ Hy/ 2mE act ] , where if is a width of the barrier, m is the 
effective mass, and E &ct is the height of the barrier, if the reaction proceeds via tunneling. 
The width of the barrier w is fixed to 1 A. The effective mass is equal to the reduced mass 
of the reactants in a two-body addition reac tion, but depends on the details of reaction for 
more complex reactions such as abstraction (IHidaka et al. 1 120091 ) . Since our model includes 
various grain-surface reactions of which the details (such as potential energy surface) are not 
yet well understood, we use the reduced mass of the reactants as the effective mass for all 
grain-surface reactions. The activation barrier of the reaction H + CO is set to 2500 K in 
our original model. Figure [6] shows reactions in our model and E act which we assume for 
each reaction in units of K; E act is determined so that the rate coefficients relative to H + 
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CO become roughly consistent with IWatanabe fc Kouchi I (120081 ) and iHidaka et al. I ( 120091 ). 
Reactions depicted by the solid arrows are included in our original model, while the dashed 
arrows are newly added reactions. In the original model, the activation barrier for D + CO is 
set to be the same as H + CO (2500 K); due to the larger effective mass, the rate coefficient 
of D + CO is lower than that of H + CO by several orders of magnitude. In the laboratory 
experiment, on the other hand, the relative rate coefficient of D + CO is about 10 % to H 
+ CO. Hence we changed the barrier for D + CO to 1400 K in the model. Similarly the 
barrier for reaction of D atom with H 2 CO (and its deuterated isotopomers) are lowered. It 
should be noted that these values of E act are not necessarily equal to the activation barrier 
obtained from the quantum chemical calculations, since the rate of the tunneling sensitively 
depends on the shape of the potential and effective mass in reality, while we assume a simple 
rectangular barrier and adopt the reduced mass as the effective mass. We also note that the 
frequency factor is often replaced by accretion or desorption rates of H and D atoms in our 
calculation, because we adopt the modified rate. 



Taquet. Ceccarelli fc Kahane I (120121 ) calculated the D/H ratios of H2CO and CH3OH 
in pseudo-time-dependent models of nn = 10 4 — 5 x 10 6 cm -3 cm -3 , including the substi- 
tution and abstraction reactions. Their grain-surface chemistry model is more sophisticated 
than ours; they consider multiple layers and cracks in the grain mantles. The rat e coef- 
ficie nts are determ i ned r elative to the H + CO rate, referring to IWatanabe et al. I (120 101 ) 



and IHidaka et al. I (120091 ) . They showed that the D/H ratios reach higher values at higher 
densities; in the model with n H = 5x 10 6 cm" 3 and T = 10 K, the CH 2 DOH/CH 3 OH ratio 
exceeds unity at t > 10 4 yr. 



In order to compare our chemical network model with the work of lTaquet. Ceccarelli fc Kahane 

( 120121 ). we calculated a pseudo-time-dependent model with % = 5 x 10 6 cm~ 3 and T = 10 



K; Figure [7] shows the temporal variation of H 2 CO, CH 3 OH and their deuterated isotopes 
in our model. The initial abundances are set by solving the network under molecular-cloud 
conditions («h = 2 x 10 4 cm -3 , T = 10 K) for 8 x 10 4 yr in panel (a) and for 2 x 10 5 yr in 
panel (b). The molecular D/H ratios become very hig h 0.1 — 1.0 at several 10 4 — 10 5 y r. Th e 



rise in the D/H ratio is slower than in the model of lTaquet. Ceccarelli &: Kahane I (120121 ). 
in which the HDCO/H 2 CO ratio exceeds 0.1 at ~ 10 3 yr and the CH 2 DOH/CH 3 OH ratio 
exceeds unity at ~ 10 4 yr. The slow rise of the D/H ratio i n our model would be mainly due 
to the initial abundances. iTaquet. Ceccarelli fc Kahane I ( 120121 ) calculate the atomic D/H 
ratio in steady state as a function of density, and adopt it as an initial abundance. The 
steady-state atomic D/H ratio is greater under higher- density conditions. Our initial abun- 
dances, on the other hand, are determined by the molecular-cloud calculation of n-^ = 2 x 10 4 
cm -3 . In our pseudo-time-dependent model of % = 5 x 10 6 cm -3 , the atomic D/H ratio 
increases in 10 4 yr to reach the steady state. Initial abundances of CO, CH3OH and H 2 CO 
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would also be relevant; the temporal variation of the molecular D/H ratio differs in panel 
(a) and (b). Although deuterated H 2 CO and CH 3 OH increase rapidly both in (a) and (b) 
at t > 10 4 yr, when the atomic D/H ratio reaches steady state, the molecular D/H ratios 
do not reach unity in pan el (b), due to the high abundance of normal isotopes in the initial 



conditions. The model of iTaquet. Ceccarelli fc Kahane I ( 120121 ). on the other hand, assumes 



that all carbon initially resides in CO. 



Having tested our chemical network against ITaquet. Ceccarelli fc Kahane I (120121 ). we 
now calculate molecular abundances and D/H ratios including the substitution and abstrac- 
tion reactions at t = —5.6 x 10 2 yr, 4.3 x 10 2 yr and 9.3 x 10 4 yr in our star-forming cores. 
The result is almost the same as in Figure 3; the D/H ratio increases at most by a factor of 
a few. There are three possible reasons for this low D/H ratio in spite of the newly-added 
reactions: 

(i) Our model starts from a dense molecular cloud core. A significant amount of CH 3 OH 
and H 2 CO is already formed during the pre-collapse phase with < 5 x 10 4 cm -3 , 
during which the deuterium enrichm ent by the new reactions are not very efficient 
( Taquet. Ceccarelli fc Kahane I I2OI2I ) . 



(ii) The period spent at high density (> 10 6 cm -3 ), in which the deuterium enrichment 
is efficient, is limited by the dynamical model. The temporal variation of density in 
the inf ailing fluid parcels that reach R = 2.5 AU at 9.3 x 10 4 yr is shown in Figure 1. 
The parcel falls onto the protostar in several 10 3 yrs after it enters the region where 
n H > 10 6 cm" 3 . 

(iii) The duration of the high-density phase is longer in fluid parcels that reach the central 
region in earlier stages (e.g. t = —5.6 x 10 2 yr and 4.3 x 10 2 yr); but the dust tempera - 



ture for such parcels is mostly very low; it could reach ~5K (e.g. lKeto fc Casellill2010l ). 
Since we assume thermal hopping, even H and D atoms cannot efficiently migrate on 
grains. 

Our model indicates that the D/H ratios of CH3OH and H2CO in the star- forming core 
depend significantly on the dynamical model, the heat balance (e.g. cosmic-ray flux) inside 
the prestellar cores, and the efficiency of grain-surface migration of H and D atoms at very 
low temperature. For example, we found that the gaseous CH 2 DOH/CH 3 OH ratio reaches 
17 % inside the sublimation radius (r < 7.5 AU) at t = 4.3 x 10 2 yr, if we allow H and 
D atoms to migrate via tunneling. Recent laboratory work excl udes the migration of H 



and D atoms via tunneling, but finds fast thermal hopping at 8 K (IWatanabe et al. Il2010l ). 



Considering the strong dependence of the hopping rate on temperature, thermal hopping of 
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H and D atoms could still be inefficient at 5 K. In such a case, the Eley-Rideal mechanism, 
which is not included in our model, could be important. We postpone further investigation 
of this issue to future work. 



5. Summary 

We have investigated molecular abundances a nd D/H ratios in a star - formi ng core by 
adopting the ID radiation hydrodynamics model of iMasunaga fc Inutsuka 



Spatial distributions of assorted molecules in the gas phase and ice mantle are reported 
for a core at a time 5.6 x 10 2 yr prior to protostar formation (at which the so-called first 
core is formed), and at 4.3 x 10 2 yr and 9.3 x 10 4 yr after the protostar is born. Significant 
depletion onto grains occurs in the early stages, at radii 10 AU < r < 10 3 AU and 100 AU 
< r < 10 3 AU, for times t = —5.6 x 10~ 2 yr and 4.3 x 10 2 yr, respectively. As the core 
gets warmer, complex species, as well as simple species such as NH 3 , sublime from the grain 
surfaces near the core center. The abundances of these complex species - both in the gas 
phase and ice mantle - increase with time, while the variation is rather small for CH3CN, 
which can be formed readily at low temperatures. Methanol, which can also be formed at 
low temperatures, decreases gradually as the core evolves. We find reasonable agreement 
with observation for our model abundances at t — 9.3 x 10 4 yr and T > 100 K, except for 
CH3OCH3, which is significantly under-abundant in our ID core model. 

Unsaturated carbon chains are formed from sublimed CH4 at T ~ 25 K (r ~ 1000 
AU) at if: = 9.3 x 10 4 yr. We sho w that the c h emist ry that forms these chains is indeed the 



WCCC mechanism developed by ISakai et al.l (120081 ) to explain the chemistry in L1527 and 



IRAS15398 - 3359. Our model also reproduces the central dip for C2H and C4H observed by 



Sakai et al.l (120101 ). In the earlier stages, however, the gas-phase formation of unsaturated 
carbon chains is much less efficient due to the higher density at the CH 4 sublimation zone. 
In order for WCCC to be efficient, the gas density at the CH 4 sublimation zone should be 
relatively low, so that C + , which is another ingredient for carbon chain formation, becomes 
abundant. 

In the cold phase (T ~ 10 K), molecular D/H ratios are enhanced by exothermic ex- 
change reactions such as H 2 D + + HD. We find very high molecular D/H ratios (close to 
unity) for gas-phase species in the CO depletion zone, because CO is the main reaction 
partner of H2D + . In the central regions of the protostellar core, however, molecular D/H 
ratios are mostly ~ 10~ 2 ; they originate in the D/H ratios formed on ice mantles, which 
are accumulated from the cloud core phase (without significant CO depletion) all the way 
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through to the protostellar phase. At t — 9.3 x 10 4 yr, a significant amount of complex or- 
ganic species and carbon chains are formed via grain-surface and gas-phase reactions at warm 
temperatures (T > 30 K for organic species and T > 25 K for carbon chains). Their D/H 
ratios are ~ 10 -2 , inheriting the D/H ratio of their precursor molecules, and vary slightly 
among species. Exothermic exchange reactions are still active at temperatures of a few tens 
of Kelvin, although the rates of the backward reactions increase with temperature. HCOOH 
has a higher D/H ratio than other complex organic species; it is formed by H 2 CO + OH, 
so that HCOOD can be formed from H 2 CO + OD where OH is deuterated by the exchange 
reaction of OH + D. This exchange reaction also significantly enhances the HDO/H2O ratio 
at T ~ 50K. 

While the observed molecular D/H ratios in protostellar cores could be mostly explained 
by the combination of D/H ratios in the hot-corino region (i.e. T > 100 K) and CO depletion 
zone, the D/ H ratio of CH3OH ar e significantly underestimated in our model . Lab oratory 



experiments (INagaoka et al.l 120051 ; IWatanabe fc Kouchi 1 120081 ; iHidaka et al. 1 120091 ) found 



substitution and abstraction reactions on grain surfaces to deuterate H 2 CO and CH 3 OH; 
we have added these reactions to our model. Although the D/H ratios of CH 3 OH and 
H 2 CO can exceed unity at t > 10 5 yr in a pseudo-time-dependent model with a high density 
(riH = 5 x 10 6 cm -3 ) and low temperature (T = 10 K), the ratios increased only by a factor 
of a few in our protostellar core model. A significant amount of CH3OH ice is formed in 
the prestellar stage, in which the density is still too low for the newly-added reactions to 
be efficient. While the high density phase lasts longer in inf ailing fluid parcels in earlier 
evolutionary stage of a core, the temperature could be as low as 5 K there. Eley-Rideal 
mechanism, rather than Langmuir-Hinshelwood mechanism assumed in our model, could 
therefore be important. Further studies are needed on this issue. 

Finally, in order to see how the molecular abundances evolve after fluid parcels land on 
the disk, we performed pseudo-time-dependent models at T = 40 K and 150 K starting from 
the molecular abundances set by the ID collapse model. At T = 40 K, a significant amount of 
HCOOCH3 ice is formed by the grain-surface reactions. C 2 Hg ice becomes abundant among 
carbon chains; it is formed via grain-surface hydrogenation of C 2 H 4 , which is produced by 
a gas-phase reaction. Since C 2 H 4 can marginally freeze onto grains at 40 K, its adsorption 
onto grains works as a sink in the gas-phase chemical reaction network. At T = 150 K, C7H4 
ice becomes the dominant carbon chain, and its D/H ratio decreases with time at t > 10 5 yr. 
Meanwhile, CH3OCH3 is efficiently formed in the high-temperature gas from the sublimed 
CH3OH. It is noteworthy that CH3OCH3 has a high D/H ratio (~ 0.040), in spite of its 
formation at high temperature, because its mother molecule, CH 3 OH, is highly deuterated 
during the cold phase. It is also interesting that CH3OCH3 is under-abundant in our core 
model compared with observation, while its higher abundance in the pseudo-time dependent 
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model is in much better agreement. The pseudo-time-dependent models also give higher 
abundances for some other complex organic species such as HCOOCH 3 compared with those 
in the ID core model. Since the size of the hot corino region (T > 100 K) is close to a typical 
centrifugal radius (~ 100 AU) of cloud cores, at least a fraction of the large organic species 
might be formed in such prolonged (i.e. longer than free-fall time scale) warm chemistry 
due to partial rotational support. Such predictions should be tested via observations at high 
spatial resolution. 
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Fig. 1. — (a — c) Radial distribution of density (%), temperature and infalling velocity at 
t = —5.6 x 10 2 yr (dotted line) , 4.3 x 10 2 yr (dashed line), and 9.3 x 10 4 yr (solid line), 
where t = is defined by the birth of a protostar. (d — e) Temporal variation of density and 
temperature in fluid parcels which reach r = 2.5 AU at t = —5.6 x 10 2 yr, 4.3 x 10 2 yr, and 
9.3 x 10 4 yr. 
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Fig. 2. — Radial distribution of simple molecules (a — b), complex organic species (c), and 
carbon chain species (d) at t = —5.6 x 10 2 yr, 4.3 x 10 2 yr, and 9.3 x 10 4 yr. Solid lines and 
dashed lines depict gas-phase species and ice mantle species, respectively while the chemical 
species are distinguished by line colors. Note that the horizontal axis starts from 2 AU rather 
than 1 AU. 
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Fig. 2.— (cont.) 
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Fig. 3. — Radial distribution of mono-deuterated species and their normal isotope counter- 
parts at t = —5.6 x 10 2 yr, 4.3 x 10 2 yr, and 9.3 x 10 4 yr. The bottom panels show the radial 
distribution of the D/H ratio for assorted molecules. The thin, thick, and bold lines depict 
the D/H ratios at t = —5.6 x 10 2 yr, 4.3 x 10 2 yr, and 9.3 x 10 4 yr, respectively. 
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Fig. 3. — cont. 
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Fig. 4. — Radial distribution of mono-, multi-deuterated species and their normal isotope 
counterparts at t = —5.6 x 10 2 yr, 4.3 x 10 2 yr, and 9.3 x 10 4 yr. 
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Fig. 5. — Pseudo-time-dependent models with nn = 1.487 x 10 8 cm -3 , T = 40 K (panels a, 
b) and T = 150 K (panels c, d). The initial abundance is set by the ID model at t — 9.3 x 10 4 
yr and r = 30 AU. 
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Fig. 6. — Reaction network of H2CO and CH3OH formation on the grain surface. Solid 
arrows refer to the reactions in our original model, while the dashed arrows are newly added 
reactions. The H or D and values represent the reactants and activation barrier heights 
expediently assumed for the tunneling of our rectangular potential. For D atom reactions 
with asterisks, we lowered the activation barrier so that their rates are about 10 % of the H 
atom reactions. 



-41 - 



X 



le-4 
le-5 
le-6 
le-7 
, le-8 
le-9 
le-10 
le-11 
le-12 
le-5 
le-6 
le-7 

© le " 8 
le-9 

le-10 

le-11 




le-12 



HDCOice , 

/ / .''' / / CHD20Hice 

/ / / / 



CH2DOHice / 

'/ / 

/ /' 
/ / 
/ / 

D2COice / / 



/ 




COice 



H2CQice 



CH30Hice 




HDCOice/' 



t / / / 

t / / / 
> / / / 

'/'/ / / CHD20Hice 



CH2DOHice / 

/ / 

D2COice/ / 
/ / 
/ / 
y i 



/ 

/ 



1000 



10000 100000 

time [yr] 



le+06 



Fig. 7. — Pseudo-time-dependent model at nn = 5 x 10 6 cm -3 and T = 10 K with the 
substitution and abstraction reactions. The initial abundance is determined by calculating 
the molecular evolution at molecular cloud conditions for 8 x 10 4 yr (a) and 2 x 10 5 yr (b). 



